MODULE SRFSN_WEBALS_MOD
CONTAINS
SUBROUTINE SRFSN_WEBALS(KIDIA,KFDIA,KLON,KLEVSN,&
 & PTMST  , PFRTI   , PCIL, LDLAND  , LLNOSNOW, &
 & PSSNM1M, PWSNM1M , PRSNM1M , PTSNM1M,&
 & PSLRFL , PSSRFLTI, PAHFSTI , PEVAPTI,&
 & PSSFC  , PSSFL   , PEVAPSNW, &
 & PTSFC  , PTSFL   ,           &
 & PTSURF , PTICE, PSNOTRS ,           &
 & PAPRS  , PWSURF  ,           &
 & PSSDP3,&
 & YDSOIL , YDCST   ,           &
 & PTSN, PGSN)
 


USE PARKIND1 , ONLY : JPIM, JPRB
USE YOMHOOK  , ONLY : LHOOK, DR_HOOK, JPHOOK

USE YOS_SOIL , ONLY : TSOIL 
USE YOS_CST  , ONLY : TCST
USE EC_LUN   , ONLY : NULERR

USE ABORT_SURF_MOD
USE YOMSURF_SSDP_MOD

! (C) Copyright 2021- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.

!**** *SRFSN_WEBALS* - Snow simplified energy balance 
!     PURPOSE.
!     --------
!     THIS ROUTINE COMPUTES ENERGY BALANCE IN THE SNOWPACK

!**   INTERFACE.
!     ----------
!          *SRFSN_WEBALS* IS CALLED FROM *SURFTSTPS*.

!     PARAMETER   DESCRIPTION                                    UNITS
!     ---------   -----------                                    -----

!     INPUT PARAMETERS (INTEGER):
!    *KIDIA*      START POINT
!    *KFDIA*      END POINT
!    *KLON*       NUMBER OF GRID POINTS PER PACKET
!    *KLEVSN*    VERTICAL SNOW LAYERS

!     INPUT PARAMETERS (REAL):
!    *PTMST*      TIME STEP                                      S

!     INPUT PARAMETERS (LOGICAL):
!    *LDLAND*     LAND/SEA MASK (TRUE/FALSE)
!    *LLNOSNOW*   NO-SNOW/SNOW MASK (TRUE IF NO-SNOW)

!     INPUT PARAMETERS AT T-1 OR CONSTANT IN TIME (REAL):
!    *PCIL*       FRACTION OF LAND GLACIER COVERED
!    *PWLM1M*     SKIN RESERVOIR WATER CONTENT                 kg/m**2
!    *PSSNM1M*    TOTAL SNOW MASS IN EACH LAYER (per unit area) kg/m**2
!    *PWSNM1M*    LIQUID WATER CONTENT IN SNOW                 kg/m**2
!    *PRSNM1M*    SNOW DENSITY in each layer                   kg/m**3
!    *PTSNM1M*    TEMPERATURE OF SNOW LAYER                    K
!    *PTSURF*     TEMPERATURE OF TOP SOIL LAYER                K
!    *PWSURF*     SOIL MOISTURE OF TOP SOIL LAYER                K
!    *PSSRFLTI*   NET SHORTWAVE RADIATION AT THE SURFACE,
!                  FOR EACH TILE                                 W/M**2
!    *PSLRFLTI*   NET LONGWAVE  RADIATION AT THE SURFACE  W/M**2
!                  For earch tile   
!    *PFRTI*      TILE FRACTIONS                                 -
!    *PAHFSTI*    TILE SENSIBLE HEAT FLUX                      W/M**2
!    *PEVAPTI*    TILE EVAPORATION
!    *PEVAPSN*    EVAPORATION FROM SNOW UNDER FOREST           KG/M2/S
!    *PSSFC*      CONVECTIVE  SNOW FLUX AT THE SURFACE         KG/M**2/S
!    *PSSFL*      LARGE SCALE SNOW FLUX AT THE SURFACE         KG/M**2/S
!    *PTSFC*      Convective Throughfall at the surface        KG/M**2/S
!    *PTSFL*      Large Scale Throughfall at the surface       KG/M**2/S
!    *PSURFCOND*  THERMAL CONDUCTIVITY OF TOP SOIL LAYER
!    *ZSNOTRS*    SOLAR RADIATION FLUX INTO THE SNOWPACK      W/m**2
!    *PAPRS*      ATMOSPHERIC PRESSURE ON BOTTOM HALF LEVEL   Pa
!    

!     OUTPUT PARAMETERS AT T+1 (UNFILTERED,REAL):
!    *PTSN*        TEMPERATURE OF SNOW LAYER                    K

!    FLUXES FROM SNOW SCHEME:
!    *PGSN*       GROUND HEAT FLUX FROM SNOW DECK TO SOIL     W/M**2   (#)

!     METHOD.
!     -------
!          

!     EXTERNALS.
!     ----------
!          NONE.

!     REFERENCE.
!     ----------
!          

!     Modifications:
!     Original   ! G. Arduini    01/09/2021
!     I. Ayan-Miguez (BSC) Sep 2023 Added PSSDP3 object for spatially distributed parameters
!     Added land-ice

IMPLICIT NONE

! Declaration of arguments 
INTEGER(KIND=JPIM), INTENT(IN)   :: KIDIA
INTEGER(KIND=JPIM), INTENT(IN)   :: KFDIA
INTEGER(KIND=JPIM), INTENT(IN)   :: KLON
INTEGER(KIND=JPIM), INTENT(IN)   :: KLEVSN
REAL(KIND=JPRB)   , INTENT(IN)   :: PTMST
REAL(KIND=JPRB),    INTENT(IN)   :: PFRTI(:,:)
REAL(KIND=JPRB),    INTENT(IN)   :: PCIL(:)

LOGICAL           , INTENT(IN)   :: LDLAND(:)
LOGICAL           , INTENT(IN)   :: LLNOSNOW(:)

REAL(KIND=JPRB)   , INTENT(IN)   :: PSSNM1M(:,:)
REAL(KIND=JPRB)   , INTENT(IN)   :: PWSNM1M(:,:)
REAL(KIND=JPRB)   , INTENT(IN)   :: PRSNM1M(:,:)
REAL(KIND=JPRB)   , INTENT(IN)   :: PTSNM1M(:,:)
REAL(KIND=JPRB)   , INTENT(IN)   :: PTSURF(:)
REAL(KIND=JPRB)   , INTENT(IN)   :: PTICE(:)
REAL(KIND=JPRB)   , INTENT(IN)   :: PSNOTRS(:,:)
REAL(KIND=JPRB)   , INTENT(IN)   :: PAPRS(:)
REAL(KIND=JPRB)   , INTENT(IN)   :: PWSURF(:)

REAL(KIND=JPRB),    INTENT(IN)   :: PSLRFL(:)
REAL(KIND=JPRB),    INTENT(IN)   :: PSSRFLTI(:,:)
REAL(KIND=JPRB),    INTENT(IN)   :: PAHFSTI(:,:)
REAL(KIND=JPRB),    INTENT(IN)   :: PEVAPTI(:,:)
REAL(KIND=JPRB),    INTENT(IN)   :: PSSFC(:)
REAL(KIND=JPRB),    INTENT(IN)   :: PSSFL(:)
REAL(KIND=JPRB),    INTENT(IN)   :: PEVAPSNW(:)
!REAL(KIND=JPRB),    INTENT(INOUT):: PTSFC(:)
!REAL(KIND=JPRB),    INTENT(INOUT):: PTSFL(:)
REAL(KIND=JPRB),    INTENT(IN):: PTSFC(:)
REAL(KIND=JPRB),    INTENT(IN):: PTSFL(:)

REAL(KIND=JPRB),    INTENT(IN)   :: PSSDP3(:,:,:)

TYPE(TSOIL)       , INTENT(IN)   :: YDSOIL
TYPE(TCST)        , INTENT(IN)   :: YDCST

REAL(KIND=JPRB)   , INTENT(OUT)  :: PTSN(:,:)
REAL(KIND=JPRB)   , INTENT(OUT)  :: PGSN(:)


! Local variables (in srfsn_driver in full non linear
REAL(KIND=JPRB) :: ZSNOWF(KLON)
REAL(KIND=JPRB) :: ZRAINF(KLON)
REAL(KIND=JPRB) :: ZEVAPSN(KLON)
REAL(KIND=JPRB) :: ZHFLUX(KLON)
REAL(KIND=JPRB) :: ZFRSN(KLON)


REAL(KIND=JPRB) :: ZTSFC(KLON)
REAL(KIND=JPRB) :: ZTSFL(KLON)

! Local variables 
REAL(KIND=JPRB) :: ZDSN(KLEVSN)    ! actual snow depth
REAL(KIND=JPRB) :: ZSNHC(KLEVSN)   ! snow heat capacity 
REAL(KIND=JPRB) :: ZICE(KLEVSN)    ! snow ice content (PSSN-PWSN)
REAL(KIND=JPRB) :: ZSNCOND(KLEVSN) ! Snow thermal conductivity 
REAL(KIND=JPRB) :: ZSNCONDH(KLEVSN+1) ! THERMAL CONDUCTIVITY IN HALF LEVEL TERM 
REAL(KIND=JPRB) :: ZTA(KLEVSN),ZTB(KLEVSN),ZTC(KLEVSN),ZTR(KLEVSN)  ! TERMS TO TRI-DIAG
REAL(KIND=JPRB) :: ZTSTAR(KLEVSN)  ! New snow temperature 
REAL(KIND=JPRB) :: ZISTAR(KLEVSN)  ! New ice content 
REAL(KIND=JPRB) :: ZWSTAR(KLEVSN)  ! New liquid water content 
REAL(KIND=JPRB) :: ZLIQF(0:KLEVSN) ! LIQUID WATER FLUX 

REAL(KIND=JPRB) :: ZGSN,ZWCAP
REAL(KIND=JPRB) :: ZQ(KLEVSN),ZMELTCHANGE(KLEVSN),ZFREZCHANGE(KLEVSN),ZTMP0(KLEVSN),&
                 & ZMELTSN(KLEVSN), ZFREZSN(KLEVSN)
REAL(KIND=JPRB) :: ZGSNRES(0:KLEVSN)

REAL(KIND=JPRB) :: ZSNVCOND
REAL(KIND=JPRB) :: ZTMP1
REAL(KIND=JPRB)  :: ZTMST,ZIHCAP,ZSOILDEPTH1,ZEPSILON
REAL(KIND=JPRB) :: ZSOILRES, ZHOICE, ZSNRES, ZWHCAP
REAL(KIND=JPRB) :: ZFF, ZWU, ZLWT,  ZINVWSAT, ZLIC, ZLAMBDASAT, ZCOND, ZKERSTEN, ZSOILCOND, ZSURFCOND
REAL(KIND=JPRB) :: ZTBOTTOM(KLON)

INTEGER(KIND=JPIM) :: JL,JK,KLACT

REAL(KIND=JPRB) :: ZTMPHC, ZTMPCOND, ZSNCONDH1, ZSNCONDH2, ZCONSLWC, ZLN10

REAL(KIND=JPHOOK) :: ZHOOK_HANDLE

!! INCLUDE FUNCTIONS

#include "fcsurf.h"

!    -----------------------------------------------------------------
IF (LHOOK) CALL DR_HOOK('SRFSN_WEBALS_MOD:SRFSN_WEBALS',0,ZHOOK_HANDLE)

!    -----------------------------------------------------------------
ASSOCIATE(RHOCI=>YDSOIL%RHOCI,RHOICE=>YDSOIL%RHOICE,&
 & RTT=>YDCST%RTT,RLMLT=>YDCST%RLMLT,&
 & RALAMSN=>YDSOIL%RALAMSN, RLAMICE=>YDSOIL%RLAMICE, &
 & RFRTINY=>YDSOIL%RFRTINY, RFRSMALL=>YDSOIL%RFRSMALL, &
 & RDSNMAX=>YDSOIL%RDSNMAX, SNHCONDAV=>YDSOIL%SNHCONDAV, &
 & SNHCONDBV=>YDSOIL%SNHCONDBV, SNHCONDCV=>YDSOIL%SNHCONDCV, &
 & SNHCONDPOV=>YDSOIL%SNHCONDPOV, &
 & RKERST1=>YDSOIL%RKERST1, RKERST2=>YDSOIL%RKERST2, RKERST3=>YDSOIL%RKERST3, &
 & RLAMBDADRY=>YDSOIL%RLAMBDADRY, RLAMBDADRYM3D=>PSSDP3(:,:,SSDP3D_ID%NRLAMBDADRYM3D), & 
 & RLAMBDAICE=>YDSOIL%RLAMBDAICE, RLAMBDAWAT=>YDSOIL%RLAMBDAWAT, &
 & RLAMSAT1=>YDSOIL%RLAMSAT1, RLAMSAT1M3D=>PSSDP3(:,:,SSDP3D_ID%NRLAMSAT1M3D), RSIMP=>YDSOIL%RSIMP, &
 & RWSAT=>YDSOIL%RWSAT, RWSATM3D=>PSSDP3(:,:,SSDP3D_ID%NRWSATM3D), &
 & RCONDSICE=>YDSOIL%RCONDSICE, &
 & RLWCSWEA=>YDSOIL%RLWCSWEA , RLWCSWEB=>YDSOIL%RLWCSWEB, RLWCSWEC=>YDSOIL%RLWCSWEC)

! RLMLT (latent heat of fusion J Kg -1)
ZTMST      = 1.0_JPRB/PTMST 
ZIHCAP     = RHOCI/RHOICE  ! Ice heat capacity (J K-1 Kg-1)
ZSOILDEPTH1= YDSOIL%RDAW(1) ! 1sth Soil layer depth
ZEPSILON   = 10._JPRB*EPSILON(ZEPSILON)
ZWHCAP     = 4180._JPRB ! J K-1 Kg-1
ZLN10      = LOG(10.0_JPRB)

! Input fields (in srfsn_driver usually)
  
DO JL=KIDIA,KFDIA
! snow fraction and logical no snow
  ZFRSN(JL)=MAX(PFRTI(JL,5)+PFRTI(JL,7),RFRTINY)

  IF (LDLAND(JL)) THEN
    ! Snowfall is all redirected here ! 
    ZSNOWF(JL) = PSSFC(JL) + PSSFL(JL)
  ELSE 
    ZSNOWF(JL) = 0.0_JPRB 
  ENDIF
! Evap. is only fractional 
  ZEVAPSN(JL) = PFRTI(JL,5)*PEVAPTI(JL,5) + PFRTI(JL,7)*PEVAPSNW(JL)

    
! Rainfall is only fractional 
  ZRAINF(JL)= ZFRSN(JL)*(PTSFC(JL)+PTSFL(JL))
  ZTSFC(JL) = (1._JPRB - ZFRSN(JL) )*PTSFC(JL)
  ZTSFL(JL) = (1._JPRB - ZFRSN(JL) )*PTSFL(JL)

  IF(.NOT. YDSOIL%LESNICE)THEN
    ZTBOTTOM(JL)=PTSURF(JL)
  ELSE
    IF (LDLAND(JL))THEN
      ZTBOTTOM(JL)=PTSURF(JL)
    ELSE
      ZTBOTTOM(JL)=PTICE(JL)
    ENDIF
  ENDIF
ENDDO


DO JL=KIDIA,KFDIA
  IF (LLNOSNOW(JL)) THEN 
    PTSN(JL,:) = RTT
    PGSN(JL)   = 0.0_JPRB
    
  ELSE
! heat flux into the snowpack (not scaled by fraction, it simplifies fraction in
! heat transfer eq.
    ZHFLUX(JL)=( PFRTI(JL,5)*(PAHFSTI(JL,5)+YDCST%RLSTT*PEVAPTI(JL,5))&
              & +PFRTI(JL,7)*(PAHFSTI(JL,7)+YDCST%RLVTT*(PEVAPTI(JL,7)-PEVAPSNW(JL))&
              & +YDCST%RLSTT*PEVAPSNW(JL))&
              & +PFRTI(JL,5)*PSSRFLTI(JL,5)&
              & +PFRTI(JL,7)*PSSRFLTI(JL,7)&
              & +(PFRTI(JL,5)+PFRTI(JL,7))*PSLRFL(JL)&
              & ) !

! Soil heat conductivity (simplified) - copied from srfts
    IF (.NOT. YDSOIL%LESNICE) THEN
      ZFF=0.0_JPRB
      ZWU=PWSURF(JL)*(1.0_JPRB-ZFF)
      ZLWT=RLAMBDAWAT**ZWU
      ZINVWSAT=1.0_JPRB/RWSATM3D(JL,1)
      ZLIC=RLAMBDAICE**(RWSATM3D(JL,1)-ZWU)
      ZLAMBDASAT=RLAMSAT1M3D(JL,1)*ZLIC*ZLWT
      IF (PWSURF(JL)*ZINVWSAT > RKERST1) THEN
        ZCOND = PWSURF(JL)*ZINVWSAT
      ELSE
        ZCOND = RKERST1
      ENDIF
      ZKERSTEN=RKERST2*LOG10(ZCOND)+RKERST3
      ZSOILCOND=RLAMBDADRYM3D(JL,1)+ZKERSTEN*(ZLAMBDASAT-RLAMBDADRYM3D(JL,1))
      ZSURFCOND= (1._JPRB - PCIL(JL))*ZSOILCOND + PCIL(JL)*RCONDSICE 
  ELSE
    IF (LDLAND(JL))THEN
      ZFF=0.0_JPRB
      ZWU=PWSURF(JL)*(1.0_JPRB-ZFF)
      ZLWT=RLAMBDAWAT**ZWU
      ZINVWSAT=1.0_JPRB/RWSATM3D(JL,1)
      ZLIC=RLAMBDAICE**(RWSATM3D(JL,1)-ZWU)
      ZLAMBDASAT=RLAMSAT1M3D(JL,1)*ZLIC*ZLWT
      IF (PWSURF(JL)*ZINVWSAT > RKERST1) THEN
        ZCOND = PWSURF(JL)*ZINVWSAT
      ELSE
        ZCOND = RKERST1
      ENDIF
      ZKERSTEN=RKERST2*LOG10(ZCOND)+RKERST3
      ZSOILCOND=RLAMBDADRYM3D(JL,1)+ZKERSTEN*(ZLAMBDASAT-RLAMBDADRYM3D(JL,1))
        ZSURFCOND= (1._JPRB - PCIL(JL))*ZSOILCOND + PCIL(JL)*RCONDSICE 
    ELSE
        ZSURFCOND=RCONDSICE
      ZSOILDEPTH1=YDSOIL%RDAI(1)
    ENDIF
   ENDIF
!! Preparation
    DO JK=1,KLEVSN
      IF (PSSNM1M(JL,JK) > ZEPSILON ) KLACT=JK
    ENDDO
    DO JK=1,KLEVSN
      ZICE(JK)    = PSSNM1M(JL,JK) - PWSNM1M(JL,JK) ! Ice (Kg m-2)
      ! Limit thermally active snow depth to 1 meter of snow (for glaciers in particular)
      ZDSN(JK) = PSSNM1M(JL,JK)/ PRSNM1M(JL,JK) ! read snow depth (m)
      IF (ZDSN(JK) < RDSNMAX) THEN
        ZDSN(JK)   = PSSNM1M(JL,JK)/ PRSNM1M(JL,JK)
      ELSE
        ZDSN(JK)   = RDSNMAX
      ENDIF
    ZTMPHC=(ZICE(JK)/PRSNM1M(JL,JK)) 
    IF (ZTMPHC < RDSNMAX) THEN
      ZSNHC(JK)  = (ZIHCAP*PRSNM1M(JL,JK) * (ZICE(JK)/PRSNM1M(JL,JK))+ ZWHCAP*PWSNM1M(JL,JK) ) * ZTMST
    ELSE
      ZSNHC(JK)  = (ZIHCAP*PRSNM1M(JL,JK)*RDSNMAX + ZWHCAP*PWSNM1M(JL,JK) ) * ZTMST
    ENDIF

    ! heat conductivity from water vapor transport into the snowpack
    ZTMPCOND=(SNHCONDAV-SNHCONDBV/(PTSNM1M(JL,JK)-SNHCONDCV)) 
    IF (ZTMPCOND > 0._JPRB) THEN
        ZSNVCOND=(SNHCONDPOV/PAPRS(JL))*(SNHCONDAV-SNHCONDBV/(PTSNM1M(JL,JK)-SNHCONDCV))
    ELSE
        ZSNVCOND=0._JPRB
    ENDIF
    !  ! snow heat conductivity (total)
      ZSNCOND(JK) = FSNTCOND(PRSNM1M(JL,JK))+ZSNVCOND  ! add the thermal cond from water vapor transfer
    ENDDO
     
    ZSNCONDH(1)=ZFRSN(JL)*ZSNCOND(1) / MAX(ZEPSILON,(0.5_JPRB*ZDSN(1))*(0.5_JPRB*ZDSN(1)))   ! ACTUALY NOT USED ! (W m-2 K-1)
    DO JK=2,KLACT
      ZSNCONDH(JK)=ZFRSN(JL)*2._JPRB*(ZDSN(JK-1)*ZSNCOND(JK-1)+ZDSN(JK)*ZSNCOND(JK))/&
                  & ((ZDSN(JK-1)+ZDSN(JK))*(ZDSN(JK-1)+ZDSN(JK)))
    ENDDO
    ! Simplified basal heat conduction
     ZSNCONDH(KLACT+1)=ZFRSN(JL)*1._JPRB*(ZDSN(KLACT)*ZSNCOND(KLACT)+ZSOILDEPTH1*ZSURFCOND)/&
                  & ((ZDSN(KLACT)+ZSOILDEPTH1)*(ZDSN(KLACT)+ZSOILDEPTH1))

    IF (KLACT > 1 ) THEN 
      JK=1
      ZTA(JK)=0._JPRB
      ZTB(JK)=ZSNHC(JK)+ZSNCONDH(JK+1)
      ZTC(JK)=-ZSNCONDH(JK+1)
      ZTR(JK)=(ZHFLUX(JL)-PSNOTRS(JL,JK))+ZSNHC(JK)*PTSNM1M(JL,JK)
! Add radiative flux to second to bottom layer:
      JK=KLACT
      ZTA(JK)=-ZSNCONDH(JK)
      ZTB(JK)=ZSNHC(JK)+ZSNCONDH(JK)+ZSNCONDH(JK+1)
      ZTC(JK)=0._JPRB
      IF(.NOT. YDSOIL%LESNICE)THEN
        ZTR(JK)=ZSNHC(JK)*PTSNM1M(JL,JK)+ZSNCONDH(JK+1)*PTSURF(JL)+PSNOTRS(JL,JK)
      ELSE
        ZTR(JK)=ZSNHC(JK)*PTSNM1M(JL,JK)+ZSNCONDH(JK+1)*ZTBOTTOM(JL)+PSNOTRS(JL,JK)
      ENDIF
      DO JK=2,KLACT-1
        ZTA(JK)=-ZSNCONDH(JK)
        ZTB(JK)=ZSNHC(JK)+ZSNCONDH(JK)+ZSNCONDH(JK+1)
        ZTC(JK)=-ZSNCONDH(JK+1)
        ZTR(JK)=ZSNHC(JK)*PTSNM1M(JL,JK) + PSNOTRS(JL,JK)
      ENDDO
    ENDIF
    
! SOLVE
    ZTSTAR(:) = PTSNM1M(JL,:)
    ZGSN=0._JPRB
    IF (KLACT == 1 ) THEN
      ! SINGLE LAYER
      ZTB(1)=ZSNHC(1)+ZSNCONDH(2)
      IF(.NOT. YDSOIL%LESNICE)THEN
        ZTR(1)=(ZHFLUX(JL)-PSNOTRS(JL,1))+ZSNHC(1)*PTSNM1M(JL,1)+ZSNCONDH(2)*PTSURF(JL)
      ELSE
        ZTR(1)=(ZHFLUX(JL)-PSNOTRS(JL,1))+ZSNHC(1)*PTSNM1M(JL,1)+ZSNCONDH(2)*ZTBOTTOM(JL)
      ENDIF
      ZTSTAR(1)=ZTR(1)/ZTB(1)
      IF(.NOT. YDSOIL%LESNICE)THEN
        ZGSN=ZSNCONDH(2)*(ZTSTAR(1)-PTSURF(JL))
      ELSE
        ZGSN=ZSNCONDH(2)*(ZTSTAR(1)-ZTBOTTOM(JL))
      ENDIF
      
    ELSE
      ! MULTI LAYER
      CALL TRISOLVERS(KLACT,ZTA,ZTB,ZTC,ZTR,ZTSTAR)
      IF(.NOT. YDSOIL%LESNICE)THEN
        ZGSN=ZSNCONDH(KLACT+1)*(ZTSTAR(KLACT)-PTSURF(JL))
      ELSE
        ZGSN=ZSNCONDH(KLACT+1)*(ZTSTAR(KLACT)-ZTBOTTOM(JL))
      ENDIF
    ENDIF 
    ! * SPECIAL CASE WHEN THERE IS MELTING ON THE 1ST LAYER
    IF ( ZTSTAR(1) > RTT ) THEN 
      IF (KLACT == 1 ) THEN
        ! SINGLE LAYER

        IF(.NOT. YDSOIL%LESNICE)THEN
          ZGSN=ZSNCONDH(2)*(RTT-PTSURF(JL))
        ELSE
          ZGSN=ZSNCONDH(2)*(RTT-ZTBOTTOM(JL))
        ENDIF
        ZTSTAR(1)=(ZHFLUX(JL)-PSNOTRS(JL,1)+ZSNHC(1)*PTSNM1M(JL,1)-ZGSN)/ZSNHC(1)
      ELSE
        ! MULTI LAYER
        ! SOLVE SYSTEM SETTING T1 == RTT 
        ZTB(1) = 1._JPRB
        ZTC(1) = 0._JPRB 
        ZTR(1) = RTT 
        ZTSTAR(:) = PTSNM1M(JL,:)  
        CALL TRISOLVERS(KLACT,ZTA,ZTB,ZTC,ZTR,ZTSTAR)
        IF(.NOT. YDSOIL%LESNICE)THEN
          ZGSN=ZSNCONDH(KLACT+1)*(ZTSTAR(KLACT)-PTSURF(JL))
        ELSE
          ZGSN=ZSNCONDH(KLACT+1)*(ZTSTAR(KLACT)-ZTBOTTOM(JL))
        ENDIF
        ! EXPLICIT 1ST LAYER TEMP CALCULATION
        ZTSTAR(1)=(ZHFLUX(JL)-PSNOTRS(JL,1)+ZSNHC(1)*PTSNM1M(JL,1)-ZSNCONDH(2)*(ZTSTAR(1)-ZTSTAR(2)))/ZSNHC(1)
      ENDIF
    ENDIF
    
!! UPDATE SOLID/LIQUID CONTENTS OF FIRST LAYER BEFORE FURTHER CALCULATIONS
    ZISTAR(:) = ZICE(:)
    ZWSTAR(:) = PWSNM1M(JL,:)
    ZLIQF(:)  = 0._JPRB
      !! INCLUDE SNOWFALL AND RAINFALL INTERCEPTION 
    ZISTAR(1) = ZICE(1)+PTMST*(ZSNOWF(JL)+ZEVAPSN(JL))
    ZLIQF(0)  = PTMST*ZRAINF(JL)
    
    !! PHASE CHANGES AND LIQUID WATER BALANCE STARTING IN LAYER 1
    ZMELTSN(1:KLEVSN) = 0._JPRB
    ZFREZSN(1:KLEVSN) = 0._JPRB
    ZGSNRES(0:KLEVSN) = 0._JPRB
    ZQ(1:KLEVSN)      = 0._JPRB
    ZTMP0(1:KLEVSN)       = 0._JPRB
    ZMELTCHANGE(1:KLEVSN) = 0._JPRB
    ZFREZCHANGE(1:KLEVSN) = 0._JPRB

    DO JK=1,1
      ! UPDATE LIQUID WATER FROM FLUX ABOVER 
      ZWSTAR(JK) = PWSNM1M(JL,JK) + ZLIQF(JK-1)
      
      ! PHASE CHANGES 
      ZTMP0(JK)   = ZSNHC(JK)*(ZTSTAR(JK)-RTT)
      !ZMELTSN = MAX(0._JPRB, MIN( ZTMP0 , RLMLT*ZTMST*ZISTAR(JK) ) )
      !ZFREZSN = MIN(0._JPRB, MAX( ZTMP0 , -RLMLT*ZTMST*ZWSTAR(JK) ) )
      ZMELTCHANGE(JK) = RLMLT*ZTMST*ZISTAR(JK) 
      ZFREZCHANGE(JK) = -RLMLT*ZTMST*ZWSTAR(JK)
      IF (ZTMP0(JK) < ZMELTCHANGE(JK) ) THEN
         ZMELTSN(JK)=ZTMP0(JK)
      ELSE
         ZMELTSN(JK)=RLMLT*ZTMST*ZISTAR(JK)
      ENDIF
      IF (ZMELTSN(JK)>0._JPRB) THEN
         ZMELTSN(JK)=ZMELTSN(JK)
      ELSE
         ZMELTSN(JK)=0._JPRB
      ENDIF
      IF (ZTMP0(JK) > ZFREZCHANGE(JK)) THEN
        ZFREZSN(JK) = ZTMP0(JK)
      ELSE
        ZFREZSN(JK) = -RLMLT*ZTMST*ZWSTAR(JK)
      ENDIF
      IF (ZFREZSN(JK) > 0._JPRB) THEN
        ZFREZSN(JK) = 0._JPRB
      ELSE 
        ZFREZSN(JK) = ZFREZSN(JK)
      ENDIF

      ZQ(JK)      = ZMELTSN(JK) + ZFREZSN(JK)
      
      ! FINAL TEMP. UPDATE
      ZTSTAR(JK) = ZTSTAR(JK) - (ZQ(JK)+ZGSNRES(JK-1))/ZSNHC(JK)
      IF ( ZTSTAR(JK) > RTT ) THEN
        ZGSNRES(JK)=ZSNHC(JK)*(RTT-ZTSTAR(JK))
        ZTSTAR(JK)=RTT
      ELSE
        ZGSNRES(JK)=0._JPRB
      ENDIF
      
      ! UPDATE SOLID / LIQUID MASS
      ZWSTAR(JK) = ZWSTAR(JK) + ZQ(JK)/RLMLT*PTMST 
      ZISTAR(JK) = ZISTAR(JK) - ZQ(JK)/RLMLT*PTMST
      
      ! LIQUD WATER BUDGET
      ZCONSLWC=MAX(0._JPRB,(RLWCSWEA-PRSNM1M(JL,JK))/RLWCSWEA)
      ZWCAP  = ZISTAR(JK)*(RLWCSWEB+(RLWCSWEC-RLWCSWEB)*ZCONSLWC)
      ZLIQF(JK) = (ZWSTAR(JK)-ZWCAP)
      IF (ZLIQF(JK)>0._JPRB) THEN
        ZLIQF(JK) = (ZWSTAR(JK)-ZWCAP)
      ELSE
        ZLIQF(JK) = 0._JPRB
      ENDIF
    ENDDO
    DO JK=KLACT+1,KLEVSN
      ZTSTAR(JK) = ZTSTAR(KLACT)
    ENDDO  
      
    !! FINAL VALUES
    DO JK=1,KLEVSN
      IF (ZTSTAR(JK)<RTT) THEN
        PTSN(JL,JK)= ZTSTAR(JK)
      ELSE
        PTSN(JL,JK)= RTT
      ENDIF
    ENDDO
    PGSN(JL) = ZGSN - ZGSNRES(KLACT) + PSNOTRS(JL,KLACT+1)

    
    IF (ANY(PTSN(JL,:)>RTT+ZEPSILON)) THEN
      write(*,*) 'tsn is above zero C'
      !*CALL ABORT_SURF('ALL SNOW MELTED')
    ENDIF
    DO JK=1,KLEVSN
      IF ((PTSN(JL,JK)<100._JPRB)) THEN
        write(NULERR,*) 'Very cold snow temperature, webals'
        write(NULERR,*) 'Tsn-1',PTSNM1M(JL,:)
        write(NULERR,*) 'Tsn',PTSN(JL,:)
        write(NULERR,*) 'SWE-1',PSSNM1M(JL,:)

        PTSN(JL,JK)=100._JPRB
        !*CALL ABORT_SURF('Very snow cold temperature')
      ENDIF
    ENDDO
      
  ENDIF 
ENDDO
                           
  
END ASSOCIATE
!    -----------------------------------------------------------------
IF (LHOOK) CALL DR_HOOK('SRFSN_WEBALS_MOD:SRFSN_WEBALS_MOD',1,ZHOOK_HANDLE)

END SUBROUTINE SRFSN_WEBALS

SUBROUTINE TRISOLVERS(KSNLAC,ZA,ZB,ZC,ZF,ZTSTAR)
USE PARKIND1, ONLY : JPIM, JPRB
IMPLICIT NONE

!! SINGLE POINT TRIDIAGONAL SOLVER ! 

!* DECLARATION OF ARGUMENTS
INTEGER(KIND=JPIM)  ,INTENT(IN)  :: KSNLAC ! NUMBER OF ACTIVE SNOW LAYERS
REAL(KIND=JPRB),INTENT(IN),DIMENSION(:)  :: ZA,ZB,ZC,ZF
REAL(KIND=JPRB),INTENT(INOUT),DIMENSION(:) :: ZTSTAR

!* LOCAL VARIABLES  
INTEGER(KIND=JPIM)  :: JK
REAL(KIND=JPRB)     :: BET
REAL(KIND=JPRB),DIMENSION(SIZE(ZB)) :: GAM

BET=ZB(1)

ZTSTAR(1)=ZF(1)/BET
DO JK=2,KSNLAC
  GAM(JK)=ZC(JK-1)/BET
  BET=ZB(JK)-ZA(JK)*GAM(JK)
  ZTSTAR(JK)=(ZF(JK)-ZA(JK)*ZTSTAR(JK-1))/BET
ENDDO
DO JK=KSNLAC-1,1,-1
  ZTSTAR(JK)=ZTSTAR(JK)-GAM(JK+1)*ZTSTAR(JK+1)
ENDDO

END SUBROUTINE TRISOLVERS


END MODULE SRFSN_WEBALS_MOD
